Adaptively Biased Molecular Dynamics for Free Energy Calculations 
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We present an Adaptively Biased Molecular Dynamics (ABMD) method for the computation of the free 
energy surface of a reaction coordinate using non-equilibrium dynamics. The ABMD method belongs 
to the general category of umbrella sampling methods with an evolving biasing potential, and is 
inspired by the metadynamics method. The ABMD method has several useful features, including a 
small number of control parameters, and an 0(t) numerical cost with molecular dynamics time t. 
The ABMD method naturally allows for extensions based on multiple walkers and replica exchange, 
where different replicas can have different temperatures and/or collective variables. This is beneficial 
not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of 
useful information that may be obtained from a given simulation. The workings of the ABMD method 
are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated 
environment. 



I. 



INTRODUCTION 



When investigating the equilibrium properties of a 
complex polyatomic system, it is customary to identify 
a suitable reaction coordinate c(ri, . . . , rjy) : R 3N i— ► Q 
that maps atomic positions , . . . , r onto the points 
of some manifold Q, and then to study its equilibrium 
probability density: 



p(0 = (*[£- tr(ri 



(angular brackets denote an ensemble average) . The den- 
sity p(£) provides information about the relative stability 
of states corresponding to different values of £ along with 
useful insights into the transitional kinetics between var- 
ious stable states. In practice, the Landau free energy^ 

/(0 - -k B Tlnp(0, 

is typically preferred over p(£), because it tends to be 
more intuitive. Either or /(£) is said to provide 
a coarse-grained description of the system — in terms 
of £ alone — with the rest of the degrees of freedom of 
the original system integrated out. Quite naturally, the 
reaction coordinate (often also referred to as collective 
variable or order parameter) is typically chosen to repre- 
sent the slowest degrees of freedom of the original system, 
although this is not formally required. 

In the past few years, several methods targeting the 
computation of /(£) using non-equilibrium dynamics 
have become popular. First methods that introduced 
a time evolving potential to bias the original potential 
energy were the the Local Elevation Method (LEM)2, by 
Huber, Torda and van Gunsteren in the MD context and 
the Wang-Landau approach in MC one 3 -. More recent ap- 
proaches include the adaptive-force bias method 4 -, and 
the non-equilibrium metadynamics^ method. These 
methods all estimate the free energy of the reaction co- 
ordinate from an "evolving" ensemble of realizations 7,8 , 
and use that estimate to bias the system dynamics, so as 
to flatten the effective free energy surface. Collectively, 



they can all be considered as umbrella sampling meth- 
ods, with an evolving potential. In the long time limit, 
the biasing force is expected to compensate for the free 
energy gradient, so that the biasing potential eventually 
reproduces the free energy surface. 

In this work, we present an Adaptively Biased Molecular 
Dynamics (ABMD) method whose implementation is par- 
ticularly efficient and suited for free energy calculations. 
The method has an 0(t) scaling with molecular dynamics 
time t and is characterized by only two control parame- 
ters. In addition, the method allows for extensions based 
on multiple walkers and replica exchange for both temper- 
ature and/or the collective variables. The ABMD method 
has been implemented in the AMBER software package 9 , 
and is to be distributed freely. 

Before discussing ABMD, it is helpful to review the 
salient features of the metadynamics (MTD) method. Es- 
sentially, the MTD method is built upon the LEM method 
by exploiting Car-Parrinello dynamics: the phase space 
of the system is extended to include additional dynami- 
cal degrees of freedom harmonically coupled to the col- 
lective variable. These additional degrees of freedom are 
assumed to have masses associated with them, and evolve 
in time according to Newton's laws. The masses are 
supposed to be large enough, so that the dynamics of 
these extra-variables is driven by the free energy gradi- 
ent. Their trajectory is then used to construct a history- 
dependent biasing potential by means of placing many 
small Gaussians along the trajectory. When combined 
with Car-Parrinello ab-initio dynamics, MTD has been 
successfully used to explore complex reaction pathways 
involving several energy barrier o 10 ! 11 ' 12 ! 13 ' 14 ! 15 ' 16 ' 17 ' 18 . 

While MTD continues to be used successfully, there are 
several known limitations associated with the initial im- 
plementation of the method, which provided the motiva- 
tion for the development of the ABMD method. First, in 
order to calculate reliable free energies with a controllable 
accuracy, long runs are needed, especially for the "cor- 
rective" follow-up at equilibrium 1 ^.. This is especially 
true for biomolecular systems, which typically are char- 
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acterized by many degrees of freedom and non-negligible 
entropy contributions to the free energies. Long runs, 
however, may be precluded by the MTD method, because 
of its unfavorable scaling with molecular dynamics MD 
time t. While one can readily speed up the original MTD 
method using such tricks as truncated Gaussians and 
kd-trees^, the bottleneck there is the explicit calcula- 
tion of the history-dependent potential. Since at every 
MD step, Gaussians from all previous time steps need to 
be added, the number of Gaussians grows linearly with 
t. The numerical cost of MTD therefore grows as 0(t 2 ) 
which, in some cases, may prove itself to be prohibitively 
expensive, especially when long runs are needed. Another 
undesirable feature is that the MTD method (at least in 
its original implementation) is characterized by a rela- 
tively large number of parameters (e.g., the masses and 
spring constants associated with the collective variable, 
the characteristics of the Gaussians to be added, multi- 
ple control parameters, etc), all of which influence the 
dynamics in an entangled and non-transparent way. A 
successful MTD run often requires a careful balancing of 
these parameters, which is especially nontrivial for multi- 
dimensional collective variables. More recent implemen- 
tations of MTD 20 have reduced the number of parameters. 
As will be discussed, the ABMD method is characterized 
by only two control parameters and scales as 0(t) with 
simulation time. 



II. THE ADAPTIVELY BIASED MOLECULAR 
DYNAMICS METHOD. 

The ABMD method is formulated in terms of the follow- 
ing set of equations: 

d 2 r d 
m *^T = Fa ~ aT U [ t \ <J (ri, ■ ■ ■ , rjv) ] , 



dr a 



dU{t\£) _ k B T 



dt 



TF 



G[£ - (? fa, . . . ,r N )], 



where the first set represents Newton's equations that 
govern ordinary MD (temperature and pressure regula- 
tion terms are not shown) augmented with the additional 
force coming from the time-dependent biasing potential 
U(t\£) (with U(t = 0|£) = 0), whose time evolution is 
given by the second equation. In the following, we re- 
fer to tf as flooding timescale and to G(£) as kernel (in 
analogy to the kernel density estimator widely used in 
statistics 21 ). The kernel is supposed to be positive defi- 
nite (G(£) > 0) and symmetric (£?(-£) = G(£)). It can 
be perceived as a smoothed Dirac delta function. For 
large enough tf and small enough width of the kernel, 
the biasing potential J7(t|£) converges towards —/(C) as 
t ->ooZ»£. 

Our numerical implementation of the ABMD method in- 
volves the following. We stick with Q = Qi x ■ ■ • x Q D 
where is either [a, b] 6 I 1 or a one-dimensional torus, 



and use cubic B-splines (or products of thereof for D > 1) 
to discretize U(t\£) in Q: 

U(t\$= ]T U m (t)B{t/AS-m), 



(2-|£l) 3 /6, 



B(0 = < £*(|£|-2)/2 + 2/3, 0<|£|<1, 



0. 



otherwise. 



We use the biweight kernel 2 - 1 for G(£): 
G(0 



48 
41 



(i-e/i) , -2<£<2 

0, otherwise 



and an Euler-like discretization scheme for the time evo- 
lution of the biasing potential: 

U m (t + At) = U m {t) + At—Gla/At; - ml , 

where a = er(ri, . . . , r^r) is at time t. Note that this time 
discretization may be readily improved. This, however, is 
not really important here, since the goal is not to recover 
the solution of the ABMD equations per se, but rather to 
flatten Z7(i|£) + /(£) in the t — > oo limit. Note also, that 
the numerical cost of evaluation of the time-dependent 
potential is constant over time, and so ABMD scales triv- 
ially as 0(t), which is computationally quite favorable. 
The storage requirements of the ABMD are also quite rea- 
sonable, especially if sparse arrays are used for U m . In 
addition, it is characterized by only two control parame- 
ters: the flooding timescale tf and the kernel width 4A£. 

ABMD admits two important extensions. The first is 
identical in spirit to the multiple walkers metadynam- 
icsZ*£^. It amounts to carrying out several different MD 
simulations biased by the same E/(£|£), which evolves via: 



«-^£G[£-. (r? , 



dt 



where a labels different MD trajectories. A second exten- 
sion is to gather several different MD trajectories, each 
bearing its own biasing potential and, if desired, its 
own distinct collective variable, into a generalized en- 
semble for "replica exchange" with modified "exchange" 
rule o 23 i 24 ' 25 . Both extensions are advantageous and lead 
to a more uniform flattening of U(t\£) + /(£) in Q. This 
enhanced convergence to /(£) is due to the improved sam- 
pling of the "evolving" canonical distribution. 

We have implemented the ABMD method in the AMBER 
package 9 , with support for both replica exchange and 
multiple- walkers. In pure "parallel tempering" replica 
exchange (same collective variable in all replicas), N r 
replicas are simulated at different temperatures T n , n — 
l,...,N r . Each replica has its own biasing potential 
U n (t\{;), n = 1, . . . , N r , that evolves according to its dy- 
namical equation. Exchanges between neighboring repli- 
cas are attempted at a prescribed rate, with an exchange 
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FIG. 1: The Ace-GGPGGG-Nme peptide in a /3-hairpin confor- 
mation (sketch). 



probability given by-2^: 
w(m\n) = 



1, A<0, 
exp(-A), A > 0, 



A 



1 



1 



+ 



1 



EI! 



e; 



(2.1) 



(2.2) 



kBT m 
1 



u m {C)-u m {C) 

U n (C)-U n (£ m ) 



where i5 p denotes the atomic potential energy. The bi- 
asing potentials are temperature-bound and converge in 
the t — > oo limit to the free energies at their respective 
temperatures. 

We have also implemented a more general replica ex- 
change scheme, where different replicas can have different 
collective variables and/or temperatures, and can experi- 
ence either an evolving or a static biasing potential (the 
latter obviously includes the case of U n (t\£) = 0). Ex- 
changes between random pairs of replicas are then tried 
at a prescribed rate. This method is simply a general- 
ization 25 of the "Hamiltonian replica exchange" method 
described in Ref[H, and reduces to it when all biasing 
potentials are static. The big advantage here is that, 
by using replicas with different collective variables, it is 
possible to obtain several one- or iwo-dimensional projec- 
tions of the free energy surface for the corresponding vari- 
ables. This is very useful because it not only increases the 
amount of information that can be gathered from a given 
simulation, but it also allows for previously obtained in- 
formation for a collective variable to be used to compute 
the free energy associated with different variables. For 
instance, suppose that in the course of a simulation it 
becomes apparent that one wishes to address additional 
questions involving different collective variables. Instead 
of starting from "scratch" , one can re-use the already ob- 
tained biasing potentials and thereby greatly accelerate 
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FIG. 2: RMS error of the free energy over 3. 3 A < R g < 6. 3 A 
at T = 300A" for the ABMD (blue) and reference MTD (red) 
simulations. 



the free energy calculation for the new variables. It is 
also worth noting that for replicas running at the same 
temperature, the exchange probability does not depend 
on the atomic potential energies (Eqs. (|2.1|) - (I2.2|) above). 
This implies that the number of replicas needed to main- 
tain acceptable exchange rates can be made independent 
of the solvent degrees of freedom, provided that one is 
interested in the properties of the solute only (so that 
the collective variables do not depend explicitly on the 
atomic coordinates of the solvent) and that the structure 
is adequately solvated. This can be exploited to sam- 
ple a solute with a minimum amount of solvent, and to 
accelerate the averaging over the solvent degrees of free- 
dom. These last two applications of the general replica 
exchange method are illustrated in the next section. 



III. 



CASE STUDY: A SHORT PEPTIDE. 



To illustrate the method, we have simulated the hy- 
drophobic Ace-GGPGGG-Nme peptide (sketched in Fig{TJ) 
in the gas phase and in a solvated environment, using cy- 
clohexane as the explicit solvent. The free energy of this 
peptide at T = 300 K in gas phase has previously been 
investigated with the MTD method 19 , and is characterized 
by a "double- well" structure (see Fig[5|), with the wells 
corresponding to the peptide in a "globular" (left mini- 
mum in the FigfS]) and a /3-hairpin (right minimum in the 
FigH]) folded conformation, respectively. While simple 
enough, the molecule possesses all the typical features of 
larger peptide systems usually studied with biomolecular 
simulations. Simulation parameters are as in a previous 
study 19 : the atoms are described by the 1999 version of 
the Cornell et al. force field 2 ^, with no cutoff for the 
non-bonded interactions. The Berendsen thermostat is 
chosen with Tt p = lfs for temperature control. The MD 
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FIG. 3: RMS error of the free energy over 3.3A < R a < 6.3A 
for ABMD simulations at T = WOK with t f = 180 ps (1), 
360 ps (2), 720 ps (3) and 1440 ps (4). 
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FIG. 4: RMS error of the free energy over 3.3A< R g <6.3A 
for multiple walkers ABMD simulations with tf — 1 ns using one 
(1), two (2), four (4) and eight (8) trajectories at T = 300K. 
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time step (At) is 1 fs for the parallel tempering simula- 
tions, and 2 fs otherwise. 

The radius of gyration of the heavy atoms was chosen 
to be the collective variable: 

R a = E — ( r « - R ^) 2 • 

a 

Here, Re = ^ a {m a /rn^) r a is the center of mass, with 
m-£ = ^2 a Tia, and the sum runs over all atoms except 
hydrogen. The initial configuration is the fully unfolded 
peptide. A reference free energy profile, whose error— 
in the region of interest is less than 0.15 kcal/mol, was 
computed for benchmarking purposes (see Appendix [A] 
for details). As a measure of the RMS free energy error, 
the following construction was used: 




where 

b 

A = ^/de(/ 1 (e)-/ 2 (0), 

a 

accounts for the arbitrary additive constants in the free 
energies fi^(0- Here a = 3.3A and b = 6.3A, which 
correspond to the physical region of interest (see Figj8]). 

Figure [2] presents the time dependence of the RMS 
free energy error for the ABMD and reference MTD rur>i2. 
Both simulations have exactly the same kernel width 
4A£ = 0.25 A, and flooding timescale Tp = 90 ps that 
corresponds to the a posteriori hills acceptance rate re- 
ported in ReflijJ. It is evident that the AMBD run is 



more accurate than the corresponding MTD run. The 
ABMD method owes its better convergence to the smoother 
time evolution of the biasing potential and accurate dis- 
cretization in Q. The amount of memory used by the 
ABMD simulation to store U m values was only ss 200 x 8 
bytes (considering double precision) for roughly 10 s tiny 
"hills" that were accumulated by the end of the run. The 
reference MTD simulation with merely 5 x 10 3 Gaussians 
required roughly 25 times more memory for the biasing 
potential (with only the positions of the Gaussians stored 
explicitly). One can expect, that ABMD will be even more 
economical when it comes to dealing with multidimen- 
sional collective variables, provided that sparse arrays 
are used for U m with only non-zero elements being stored 
explicitly. 

Although an a priori error estimate for this type of 
non-equilibrium simulation is really not feasible, it is ex- 
pected that the error should decrease for increased Tp. 
This point is illustrated in Fig|3j which shows the error 
for increasing values of Tp . 

In order to decrease the simulation time required for 
accurate free energy estimates even further, the multiple- 
walker variation of ABMD proves to be useful. For a mod- 
erate number of walkers, the speedup is nearly linear, 
with an additional increase in accuracy coming from the 
better sampling of the "evolving" canonical distribution 
(see FigHJ. Parallel tempering improves both the speed 
and the accuracy even more. To this end, we first ran 
ABMD with Tp = 90 ps using 2, 4, 6 and 8 replicas at 
T = 300 K, 331 K, 365 K, 403 K, 445 K, 492 K, 543 K and 
600 K (during equilibrium MD runs, the peptide configu- 
ration jumps between the two minima on a picosecond 
timescale at T = 600K). In all cases, the Erms was 
found to be ~ 1 kcal/mol, or less as t — ► oo (data not 
shown). Again, the improvement in accuracy stems from 
the better sampling of the "evolving" canonical distribu- 
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FIG. 5: RMS error of the free energy over 3.3A < R g < 6.3A 
at T = 300A' for parallel tempering ABMD simulations using 
eight replicas running at T = 300 K, 331 K, 365 K, 403 K, 
445 K, 492 K, 543 K and 600 K with t f = 90 ps (1), 45 ps (2), 
22.5 ps (3) and 11.25 ps (4). 
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^oh kcal/mol 

FIG. 6: Free energy map for Ace-GGPGGG-Nme peptide in the 
gas phase as a function of the collective variables R g and iVcm • 
(see Refl27l for details regarding the algorithm used to make 
this plot). 



tion. Then, we ran 8 replicas with smaller values Tp and 
were surprised that the accuracy does not degrade, even 
for Tp = 11.2 ps as shown in Fig[5] 

Finally, we turn to aspects related to the general 
replica exchange method and illustrate its potential. As 
already noted, by using replicas with different collective 
variables and swapping these at prescribed rates, it is 
possible to obtain projections of the free energy surface 
for the corresponding variables. It is also possible to 
use previously obtained information with respect to one 
collective variable to compute the free energy associated 
with a different variable. For example, suppose that in- 
stead of the one- dimensional free energy profile as as a 
function of R g already discussed, one realizes that what 
is actually needed is a two-dimensional profile that in- 
cludes information with respect to the number of O-H 
bonds along the backbone. The two-dimensional free en- 
ergy map is computationally quite expensive, but the 
calculation can be greatly accelerated with the help of 
the general replica exchange method. We therefore sim- 
ulated 8 + 1 = 9 replicas. The eight replicas were run 
at the previously stated temperatures, with each replica 
biased by a static (not evolving) biasing potential corre- 
sponding to the negated free energy associated with the 
radius of gyration R g at the corresponding temperature, 
as shown in Fig(5) The additional ninth replica was run 
at T = 300 K with ABMD flooding in the two collective 
variables, i.e., R g and the number of O-H bonds along 
the backbone as given by 

o,h 1 - { roR/ro) 
where toh is the distance between a pair of hydrogen 



and oxygen atoms and tq — 2. 5 A. The sum runs over 
the unique O-H pairs (i.e., each O-H pair is counted only 
once), with O and H separated by one or more amino- 
bases along the backbone (27 pairs in total). In other 
words, we "re- use" the previously computed free energies 
for R g to get the free energy in the two-dimensional space 
(R g , Nob), the eight first replicas serve as a "sampling 
enhancement device" for the ninth replica. The calcula- 
tion is carried out in two stages: a "coarse" stage (15 ns 
with t f = 10 ps and 4:AR g = 0.25A, 4A7V OH = 0.5) 
followed by a "fine" stage (50 ns with Tp = 100 ps and 
4AR g = O.lA, 4A7V h = 0.25). In both runs exchanges 
between four randomly chosen pairs of replicas were at- 
tempted every 100 fs. The final free energy map is shown 
in FiglH] It is clear that this two-dimensional free en- 
ergy landscape conveys additional information not con- 
tained in the one-dimensional free energy plots already 
discussed. In particular, it allows for a better character- 
ization of the "globular" states of the Ace-GGPGGG-Nme 
: specifically, it is apparent from the FigJB] that there 
are at least two such states with different values of TVoh 
(both correspond to the left minimum in FiglHJ) - Of 
course, one could have re-used the information in the 
one-dimensional R g profiles to include other collective 
variables, in addition to -/Voh- 

The general replica exchange ABMD may also be advan- 
tageous for explicit solvent simulations, which are often 
notoriously lengthy. Specifically, if one is interested in 
the solute and the collective variables do not depend on 
the solvent degrees of freedom, then the number of repli- 
cas required to maintain an adequate exchange rate de- 
pends only very weakly on the amount of solvent (which 
must of course be sufficient as to adequately solvate the 
structure), provided that all the replicas are simulated 
at the same temperature. This is because the exchange 



6 



probability does not explicitly depend on the potential 
energy difference when the temperature of the replicas is 
the same. While not every choice of collective variables 
for different replicas will lead to decent exchange rates, 
one can nevertheless take advantage of this property and 
use general replica exchange to enhance the sampling in 
a solvated environment. 

In order to demonstrate the method in this regime, we 
simulated the Ace-GGPGGG-Nme peptide at T = 300 K 
solvated by 171 cyclohexane (CqRu) molecules (the to- 
tal number of atoms was 3139) under periodic boundary 
conditions using the General^ AMBER Force-Field (GAFF) 
for the solvent. We used a truncated octahedron cell of 
fixed size (constant volume) that corresponds to the equi- 
librium density at T = 300 K (the equilibrium density 
value was obtained from a 10 ns simulation under con- 
stant pressure at T = 300 K). The Particle-Mesh Ewald 
(PME 29 ) method was used for the electrostatic forces, with 
a 36 x 36 x 36 FFT grid and an 8 A cutoff for the direct 
sum (same cutoff was used for van der Waals interac- 
tions). First, we ran 10 replicas in the "flooding" mode 
(i.e., under evolving biasing potentials) using as collec- 
tive variables the distances rcc between the backbone 
carbons separated by at least 2 amino-acids (there are 
10 such distances for Ace-GGPGGG-Nme ). We ran for 5 ns 
with Tp = 30 ps and 4Arcc = 1 A, and then for another 
5 ns with Tp = 150 ps and 4Arcc = 0.5A, attempting ex- 
changes between five randomly selected pairs every 0.5 ps. 
As expected, at the start of the simulation, the exchange 
rate was nearly 100% decreasing later as the biasing po- 
tentials were built up. By the end of the simulation, 
when all possible values of the distances had been cov- 
ered, the exchange rate was very disparate between dif- 
ferent pairs of replicas. However, for every replica there 
was at least one other replica such that the exchange rate 
between them was reasonable (i.e., the whole simulation 
did not degenerate into ten different non-interacting tra- 
jectories). We then set Tp — oo in these 10 replicas, and 
added an eleventh replica whose collective variable was 
chosen as the radius of gyration of the heavy atoms. In 
other words, as before, we use the ten replicas as a "sam- 
pling enhancement device" for the last one. We then ran 
a two-stage flooding scheme: a 5 ns coarser stage, with 
Tp = 25 ps and 4Ai? g = 0.25A; followed by a 10 ns finer 
stage, with Tp = 180 ps and 4Ai? g = 0.2A. As before, 
the exchange attempts between 5 randomly selected pairs 
of replicas were performed every 0.5 ps. The "raw" ABMD- 
computed free energy associated with R g after that stage 
is shown in FigfTJ In a next step, we ran 64 biased sim- 
ulations (each comprising of 11 replicas) for 7 ns (first 
2 ns for equilibration followed by 5 ns of "production" 
runs) starting from different initial configurations. We 
set Tp = oo in all replicas (static biasing potentials) and 
recorded the values of R g in the eleventh replica every 
10 picoseconds. We then used the log-spline algorithm of 
StoneSi at. al. to estimate the (biased) log-density of the 
R g values at equilibrium. This led us to the final shape of 
the free energy curve shown in Fig[7] Compared to the 
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FIG. 7: Free energy for Ace-GGPGGG-Nme peptide solvated 
in cyclohexane at T = 300 K as obtained via: coarse 
(non-equilibrium) replica-exchange ABMD (1); finer (non- 
equilibrium) replica-exchange ABMD (2); including the correc- 
tion coming from equilibrium biased replica exchange (3). 



gas phase, the folded /3-turn in the cyclohexane solvated 
peptide is clearly favored over the globular structure. 



IV. CONCLUSIONS AND OUTLOOK. 

In summary, we have presented an ABMD method that 
computes the free energy surface of a reaction coordinate 
using non-equilibrium dynamics. The method belongs 
to the general category of umbrella sampling methods 
with an evolving potential, and is characterized by only 
two control parameters (the flooding timescale and the 
kernel width) and a favorable 0(t) scaling with molec- 
ular dynamics time t. This scaling can be very impor- 
tant for large-scale classical MD biomolecular simulations 
when long simulation times are required (see, for example 
Refl3ll. and references therein). 

ABMD has also been extended to support multiple walk- 
ers and replica exchange. Both variations improve speed 
and accuracy of the method due to the better sampling 
of the "evolving" canonical distribution. The replica ex- 
change ABMD has been generalized to include different 
temperatures and/or collective variables, that move un- 
der either an evolving or a static biasing potential. Aside 
from enhancing the sampling, this swapping of replicas 
has several important practical advantages. Most impor- 
tantly, it enables one to obtain projections of the free 
energy surface for any number of collective variables one 
might wish to investigate. In addition, one can re-use 
previously obtained results in order to enhance the sam- 
pling of new collective variables. It is also possible to 
exploit the fact that exchange rates at the same temper- 
ature are independent of the potential energy to enhance 
sampling of a solute in a minimum amount of solvent (for 



7 




, i i i i i i i , i , i i i i 

3.25 3.5 3.75 4 4.25 4.5 4.75 5 5.25 

R a (A) 

FIG. 8: The accurate free energies of Ace-GGPGGG-Nme pep- 
tide in gas phase as function of R g at T = 300 K, 331 K, 365 K, 
403 K, 445 K, 492 K, 543 K and 600 K (from bottom to top). 

collective variables independent of solvent atom coordi- 
nates). We have implemented the ABMD method in the 
AMBER package^, and plan to distribute it freely. Here, we 
have demonstrated the workings of the ABMD method with 
a study of the folding of the Ace-GGPGGG-Nme peptide, 
The application of ABMD to more complicated biomolec- 
ular systems is reserved for future publications. 
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APPENDIX A: REFERENCE FREE ENERGY 
CURVE. 

Here, we provide simulation details with regards to the 
reference free energy curve. We first ran short (5 ns, 
eight walkers with Tp = 60ps and 4A£ = 0.2A) multiple- 
walkers ABMD at T = 600 K to reconstruct the global well. 
This was followed by parallel tempering ABMD runs, using 
the biasing potential obtained from the multiple-walkers 
simulation as the zero-time value for the biasing poten- 
tials at different temperatures. We used eight replicas 
at T = 300 K, 331 K, 365 K, 403 K, 445 K, 492 K, 543 K 
and 600 K and attempted exchanges every 100 MD steps 
(0.1 ps). The simulation started with tf — 60ps and 
4A£ = 0.2A, and ran for 10 5 exchanges. We then set 
Tp to 600ps, 4A£ to O.lA and ran for 5 x 10 5 more ex- 
changes. Finally, this was followed up with 1 x 10 6 more 
exchanges with Tp = 6ns and 4A£ = O.lA. 

We then ran a very long (3 x 10 7 exchanges, 0.1 ps 
between exchanges) biased parallel tempering simula- 
tion in the spirit of Rcf 19, in order to get an a pos- 
teriori error estimate. From the resulting histogram it 
follows that the error does not exceed « 0.15kcal/mol 
for 3.3A< R g <6.3A. The RMS error is probably much 
smaller, since 0.15 corresponds to the absolute non- 
uniformity of the histogram, i.e., the maximum error, 
over 3.3A< R g <6.3A. The accurate free energy curves 
as a function of temperature are shown in FigJSj 
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